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We study neutrinoless double beta decay of several isotopes with state-of-the-art beyond self- 
consistent mean field methods to compute the nuclear matrix elements (NME). Generating coor- 
dinate method with particle number and angular momentum projection (GCM+PNAMP) is used 
for finding mother and granddaughter states and evaluating transition operators between different 
nuclei. We analyze explicitly the role of the deformation, pairing and configuration mixing in the 
evaluation of the NME. 

Double beta decay is a process where an even-even isotope decays into a nucleus with two less (more) 
neutrons (protons) with the emission of either two electrons and two neutrinos (2z//3/3) or only two 
electrons (Ou/3/3). This process becomes the only decaying mode for those nuclei where the single beta 
decay to the odd-odd neighbor is energetically forbidden [lj. While the first mode has been already 
observed for several nuclei with half-lives of ~ 10 19-21 Vj there is only one controversial claim of detection 
of neutrinoless double beta decay [2j. This mode is very relevant because can only occur, beyond the 
standard model, if neutrinos are Majorana particles. Moreover, the inverse of the half-life of this process 
is directly related to the absolute mass of the neutrinos [I]: 



where (mpp) is the effective Majorana neutrino mass, m e is the electron mass, G i is a kinematical space 
factor and finally, M 0u is the nuclear matrix element (NME). Due to the relevance of this process in 
particle and nuclear physics, there are several experiments devoted to search for Ois/3/3 decay of different 
possible emitters [3]. From the theoretical point of view, the most important quantity to determine 
is the NME. Several nuclear structure methods have been applied to compute Ou/3/3 NME being the 
Interacting Shell Model (ISM) [U [5] and proton-neutron Quasiparticle Random Phase Approximation 
(QRPA) in different versions [HI El E] the most used ones. Additionally, calculations performed with 
two other methods - angular momentum Projected Hartree Fock Bogoliubov (PHFB) with a pairing 
plus quadrupole schematic interaction jTQl [TT] and Interacting Boson Model (IBM) [12J- have been 
reported. Recently, we have proposed to study M 0l/ matrix elements with energy density functional 
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methods including beyond mean field effects [T3] . In this contribution, we present the calculation of 
NME, total Gamow- Teller (GT) strengths £+(-) and Ikeda sum rule within this framework. We also 
show results for most of the isotopes that are considered as best candidates for detecting Ov/3/3 decay 
giving a detailed analysis for 48 Ca, 76 Ge and 150 Nd of interest for CANDLES [14J, GERDA [15j and 
SNO+ [16J experiments respectively. 

Energy density functional (EDF) methods have been extensively applied for describing properly 
many properties along the whole nuclear chart (see Ref. [T7j for a review). In this work we use an EDF 
based on the Gogny D1S interaction [18] where particle number and rotational symmetries are restored 
and shape mixing (axial quadrupole degree of freedom) is taken into account within the framework of 
the Generating Coordinate Method (GCM+PNAMP functional). Contrary to other methods, single 
particle energies and residual interactions come from the same functional self-consistently. We assume 
the closure approximation for calculating NME's because at present it is not possible to compute the 
intermediate odd-odd nucleus within this framework. Then, M 0l/ of Eq. [l]can be evaluated as pQ: 

M°" = (0-f\M° u \0f) ; M°" = - (J^j M / + M% T - M*" (2) 

where |0^) are the initial and final ground states, gy — I^Qa — 1-25 are the vector and axial vector 

coupling constants and Mf/gt/t are the Fermi, Gamow- Teller and Tensor two-body operators. In this 
work, tensor term is neglected [5j [8] and F and GT operators can be written as: 

m% = v F m*\ m°j t = v GT (l) . «? (2 y^ 2) (3) 

Here t- is the isospin ladder operator that changes neutrons to protons, a are the Pauli matrices acting 
on the spin part and (ti^Vf/gt^i^) = Vf/gt(\^i — ^|) are local potentials that depend on the relative 
coordinate of the nucleons involved in the decay, r = \r\ — f^|. The functions Vf/gt{^) are the so-called 
neutrino potentials including high order currents [6] , nucleon finite size corrections [6j and radial short 
range correlations treated within the Unitary Correlator Method [T91I8] (see Ref. [5j for further details). 
We now describe the method for calculating the ground states |0^) and observables using these wave 
functions. For the sake of simplicity, we formulate the theoretical framework in terms of operators and 
wave functions although a derivation in terms of densities can be found in Refs. [201 EI]- The starting 
point is the GCM wave function: 

|o + ) = E^ /=0 ^ z l^) (4) 

p 

where P I=0 = Jq /2 e~ lJ y b sin b db, P N = ± ^ e^-^dip, P z = ± e^-^dip are the corresponding 
angular momentum (/ = 0) -for axial symmetric wave functions- and particle number projectors [22] . 
The intrinsic axial symmetric Hartree-Fock-Bogoliubov (HFB) wave functions are solutions to the 
variation after particle number projection equations constrained to a given value of the quadrupole 
deformation, f3 [23j [24]. These intrinsic wave functions are vacuum for the quasiparticle operators a 7 
that are defined by the HFB transformation [22] : 

a\ = J2 ^t4 + V57CS (5) 
s 

with U,V being the HFB variational coefficients and c\(cs) creation (annihilation) operators in the 
arbitrary single particle basis in which we expand the many body intrinsic wave functions. In our 
case, this basis corresponds to a spherical harmonic oscillator potential solved in cartesian coordinates 
that includes eleven major oscillator shells. Now, the coefficients gp are found by solving the Hill- 
Wheeler-Griffin (HWG) equation [22|, [20j [21] . First, for each nucleus we transform the non-orthogonal 
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set of wave functions {P I=0 P N P z \$ p)} into an orthonormal one {|A) = Ep ^=P I=0 P N P z \$p)} by 
diagonalizing the norm overlap matrix, ^2p/($p\P I=0 P N P z \$p')u/ i . i pf — n A u A,p- In this basis, the HWG 
equation reads: YjK' £ JSA , G% — E a G\, where £\\f are the so-called energy kernel j2H [20l [21] . Finally, 
the coefficients for the lowest eigenvalue are used to compute both the so-called collective wave functions 
P(P) = G\u^p - probability distribution for the state to have a given deformation - and all observ- 
ables. For scalar -under rotation- operators O, matrix elements between |0^) states can be calculated 
with the expression: 



0=EE ( U ^)G%^, t \P N fP z fOP^P^P z ^, i )Gl i 

A t A f Pi/3, \V nA f ) 



u K,Pi 



(6) 



In this work we will focus our interest on calculating Ov/3/3 NME's and also ground state energies, charge 
radii and total Gamow- Teller strengths <S+/-. We will compare these quantities with the available 
experimental data to check the reliability of the method. In particular, the dependence on deformation 
of the OvfiP NME's is evaluated with the kernel of Eq. [6] for the specific case of 6 = M® v (f = F/GT): 

(<b 0f \P N fP z fM9 v P I=o P N *P z *\<S>a.) 

f y]{$f} t \ P I= °P N f P Z f | <$>p f )4Wp> I P'=°P"< P* | ) 

Finally, we evaluate the total GT strength calculating the matrix elements of the operator GT = &t+/- 
between the ground state |0 + ) of the even-even initial nucleus and the 1+ states in the odd-odd final 
nucleus: 

s +/ - = E I<i^gtIo + }| 2 = (o + |(^ T )t • d^ T |o + ) (8) 

1+ 

In the r.h.s of Eq. [8] we have taken into account the completeness relation of the final states - 
1 1™) Ural = 1" th a ^ leads to a two-body form of the corresponding operators: 

S + = 3Z- 6 2h ;S. = 3N- 6 2h ; 6 2b = £ (a) a< (a) jS ft a b s a\a<; (9) 

where (<?) a £ are the matrix elements of the Pauli operators in the harmonic oscillator single particle 

basis and b a (a£) annihilates (creates) a neutron (proton) in the orbit a (£). From the above expression 
we can also see that Ikeda's sum rule (S- — S+ = 3(N — Z)) is fulfilled if particle number conserving 
wave functions are considered. 

Before evaluating the Ovfifi NME's within the framework described above, we check the reliability 
of the method comparing theoretical and experimental results for masses, charge radii and total GT 
strengths. We observe a very nice agreement between the computed values and the experimental ones. 
For example, relative errors for masses and radii are less than 1.5% in all cases. It is important to note 
that total GT strengths have been quenched by a factor (0.74) 2 as it is usually done in ISM [25] 
and QRPA [26J calculations. We now study the dependence of the Ovfifi NME's on deformation and 
pairing correlations of the isotopes involved in the process. In Figure 

|5d)- (f) we plot the magnitude of 

the GT contribution as a function of the quadrupole deformation of mother and granddaughter nuclei 
for A = 48, 76 and 150 decays (Eq. [7]) -Fermi contributions have a similar form and are not shown. 
We observe that the intensity is distributed rather symmetrically along the fa = /3f direction having 
its larger values in regions close to the spherical points. This region is wider for the A = 48 case than 
for A = 76 and 150. Apart from this area, the NME is strongly suppressed whenever the difference in 
deformation of the two nuclei is large. This has been also noticed in ISM [?] and PHFB calculations 
[TOl [TT] for 0isf3{3 and QRPA calculations [26] for 2vf3fa Nevertheless, we also obtain maxima and 
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Isotope 


BE m (MeV) 


BE ex P (MeV) [27J 


R th (fm) 


R exp ( fm ) ^8J 


Gjtheo 


gexp 


4 «Ca 


420.623 


415.991 


3.465 


3.473 


13.55 


(14.4±2.2 [50]) 


48 Xi 


423.597 


418.699 


3.557 


3.591 


1.99 


(1.9 ± 0.5 ISU]) 


76 Ge 


664.204 


661.598 


4.024 


4.081 


20.97 


(19.89 |20]) 


76 Se 


664.949 


662.072 


4.074 


4.139 


1.49 


(1.45 ± 0.07 [ST]) 


82 Se 


716.794 


712.842 


4.100 


4.139 


23.56 


(21.91 [20]) 


82 Kr 


717.859 


714.273 


4.130 


4.192 


1.24 




96 Zr 


829.432 


828.995 


4.298 


4.349 


27.63 




96 Mo 


833.793 


830.778 


4.319 


4.384 


2.56 


(0.29 ± 0.08 [33]) 


100 Mo 


861.526 


860.457 


4.372 


4.445 


27.87 


(26.69 ESI) 


100 Ru 


864.875 


861.927 


4.388 


4.453 


2.48 




116 Cd 


988.469 


987.440 


4.556 


4.628 


34.30 


(32.70 [20]) 


116 Sn 


991.079 


988.684 


4.567 


4.626 


2.61 


(i-09+°i 3 [31) 


124 Sn 


1051.668 


1049.96 


4.622 


4.675 


40.65 




i24 Te 


1051.562 


1050.69 


4.664 


4.717 


1.63 




i28 Xe 


1082.257 


1081.44 


4.686 


4.735 


40.48 


(40.08 [20]) 


128 Xe 


1080.996 


1080.74 


4.723 


4.775 


1.45 




i30 Xe 


1096.627 


1095.94 


4.695 


4.742 


43.57 


(45.90 M) 


130 Xe 


1097.245 


1096.91 


4.732 


4.783 


1.19 




136 Xe 


1143.333 


1141.88 


4.756 


4.799 


46.71 




136 Ba 


1143.202 


1142.77 


4.786 


4.832 


0.96 




i50 Nd 


1234.512 


1237.45 


5.034 


5.041 


50.32 




150 Sm 


1235.936 


1239.25 


5.041 


5.040 


1.45 





Table 1: Masses, rms charge radii and total Gamow- Teller strengths £-(+) for mother (granddaughter) 
calculated with Gogny D1S GCM+PNAMP functional compared to experimental values. Theoretical 
values for are quenched by a factor (0.74) 2 . 



minima along the diagonal Pi = f3f finding a non-trivial structure. Concerning the magnitude of the 
GT NME's, we obtain less intensity in A = 48 than in the rest, having maximum values of 3.3, 5.7 
and 5.8 for A = 48, 76 and 150 respectively. In order to shed light on this structure we represent 
in Fig. [T](g)-(i) the results for a modified GT operator where all the involved spatial dependence of 
the neutrino potentials are substituted by constant potentials -Vqt — Vgt,o = —1.5,-2.3,-2.2 for 
A = 48, 76, 150 respectively. With this assumption, GT 0is{3{3 matrix elements reduce to the 2isf3f3 ones 

evaluated in the closure approximation because oc (a -a 

^«£(2) ^ gee Eq> ui We obgerve that 

the structure as a function of the deformation is reproduced with this simplified operator indicating that 
Mqt(A, Pf) matrix elements are almost proportional to M c 2 1 zy (/5^, f3f). In addition, we can also relate the 
structure of the operator given in Fig. [j](d)-(f) with the amount of pairing correlations in the mother 
and granddaughter nuclei separately (pairing energy — E pp [22] ) represented in Fig. [l](j)-(l). We can 
observe a direct correlation between maxima and minima found both in the pairing energy and in the 
intensity of the GT contribution. Furthermore, the value of the strength is bigger with larger pairing 
correlations, in agreement with ISM calculations where zero seniority calculations give the largest values 
of the NME's [4]. 

Finally, we have to take into account the effect of configuration mixing in the calculations. In 
Fig. [j](a)-(c) we show the collective wave functions for initial and final states of A = 48,76,150 
decays respectively. Here, we observe for 48 Ca a rather constant distribution of probability between 
P & [—0.2, +0.2] with the maximum at the spherical point; for 48 Ti we find two maxima in /3 ~ —0.2 
and P ~ +0.3 and a minimum in P = giving in average a very slightly prolate deformed state; this 
is also found for 76 Ge with maxima at P ~ ±0.20 while 76 Se is mostly oblate deformed with a peak at 
P ~ —0.25. Finally, both 150 Nd and 150 Sm are well prolate deformed but their wave functions peak at 
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Figure 1: (a)-(c) Collective wave functions, GT intensity with, (d)-(f) full and, (g)-(i) constant spatial 
dependence and (j)-(l) pairing energies for (left) A = 48, (middle) A = 76 and (right) A = 150 decays. 
Shaded areas corresponds to regions explored by the collective wave functions. 

different deformations (f3 w +0.40 and f3 w +0.25, respectively). According to Eq. [6j the final results 
depend on the convolution of the collective wave functions with the matrix elements as a function 

of deformation. In Fig. [l](d)-(f ) we show schematically -shaded circles- the areas of the GT intensity 
explored by the collective wave functions. We observe, on the one hand, that configuration mixing is 
very important in the final result because several shapes can contribute to the value of NME, especially 
in A = 48 and 76. On the other hand, we see that the regions with largest values of the GT intensity 
are excluded by the collective wave functions. For example, calculations assuming spherical symmetry 
give systematically larger NME -except for A = 96- as we show in Figure [2j 

To summarize, we have presented a method for calculating nuclear matrix elements based on 

Gogny D1S Energy Density Functional including beyond mean field effects such as symmetry restoration 
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Figure 2: 0i/f3(3 matrix elements evaluated with full configuration mixing (circles) and assuming spherical 
symmetry (squares). 



and shape mixing. We have validated our method comparing theoretical and experimental ground state 
properties. Then, we have studied in detail the dependence on deformation of NME in A = 48, 76, 150 
decays relating the structure obtained with the pairing correlations in the mother and granddaughter 
nuclei. We have also found a close connection between the strength of the NME's and the results 
obtained assuming a constant dependence in the spatial part of the transition operators. Finally, we 
have pointed out the relevance of having configuration mixing in the calculations. 
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2008 and FPA2009-13377-C02-01 (MICINN). GMP is partly supported by the DFG through contract 
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